## AR(3) model as in Newell-Pizer, but with exponential smoothing/learning about r*, the long-run mean

rm(list=ls())
set.seed(616)
library(dplyr)
library(dynlm)
library(sandwich)
graphics.off()
source("R/data_fns.R")
source("R/sdr_fns.R")

##################################################
## CHOOSE HERE

## spec <- "y1"
## data <- loadShortRate() %>% rename(r = rr)

spec <- "y10"
data <- loadLongRate() %>% rename(r = rr)

##################################################

ewma <- function(x, alpha) {
    ## calculate exponentially-weighted moving average
    tau <- numeric(length(x))
    tau[1] <- x[1]
    for (t in 2:length(x))
        tau[t] <- tau[t-1] + (1-alpha)*(x[t] - tau[t-1])
    tau
}

year1 <- 1990
year2 <- 2019
N <- 400                # maximum maturity
mats <- 1:N             # yield maturities in years (for plotting)

## relevant subsamples
d1 <- data %>% filter(year <= year1)
d2 <- data %>% filter(year <= year2)

alpha <- 0.98 # smoothing parameter
data$rstar <- ewma(data$r, alpha)
cat("EWMA rstar", year1, ":", data %>% filter(year==year1) %>% pull(rstar), "\n")
cat("EWMA rstar", year2, ":", data %>% filter(year==year2) %>% pull(rstar), "\n")

## (1) take out mean - EWMA
mm <- lm(r ~ 1, weights=alpha^(-year), data=d1)
cat("Mean", year1, ":", (etabar1 <- mm$coef), "\n")
cat("SE mean:", (etase1 <- sqrt(sandwich::NeweyWest(mm, prewhite=FALSE, lag=4))[1]), "\n")
mm <- lm(r ~ 1, weights=alpha^(-year), data=d2)
cat("Mean", year2, ":", (etabar2 <- mm$coef), "\n")
cat("SE mean:", (etase2 <- sqrt(sandwich::NeweyWest(mm, prewhite=FALSE, lag=4))[1]), "\n")
d1 <- mutate(d1, rtilde = r - etabar1)
d2 <- mutate(d2, rtilde = r - etabar2)

## (2) estimate OLS on residuals
## (a) first subsample
mod1 <- dynlm(rtilde ~ L(rtilde,1) + L(rtilde,2) + L(rtilde,3) - 1, data=ts(d1), weights=tail(alpha^(year1-d1$year),-3))
print(summary(mod1))
cat("# Sum of AR coeffs pre-1990:", sum(mod1$coef), "\n")
rhobar1 <- mod1$coef
rhose1 <- sqrt(diag(vcovHC(mod1, type="HC0")))
sig1 <- summary(mod1)$sigma
## (b) second subsample
mod2 <- dynlm(rtilde ~ L(rtilde,1) + L(rtilde,2) + L(rtilde,3) - 1, data=ts(d2), weights=tail(alpha^(year2-d2$year),-3))
print(summary(mod2))
cat("# Sum of AR coeffs full sample:", sum(mod2$coef), "\n")
rhobar2 <- mod2$coef
rhose2 <- sqrt(diag(vcovHC(mod2, type="HC0")))
sig2 <- summary(mod2)$sigma

## (3) calculate term structures and save results
P1 <- sim_ar_sdr(etabar1, etase1, rhobar=rhobar1, rhose=rhose1, sig1)
P2 <- sim_ar_sdr(etabar2, etase2, rhobar=rhobar2, rhose=rhose2, sig2)
y1 <- -100/(1:N)*log(P1)
y2 <- -100/(1:N)*log(P2)
rstar_old <- etabar1
rstar_new <- etabar2
save(file=paste0("results/sdr_ar3_expon_", spec, ".RData"), rstar_old, rstar_new, P1, y1, P2, y2)

## plot term structures
dev.new()
plot(1:N, y1, type="l", ylim=range(0, y1, y2, 4), lwd=2, ylab="Discount rate (percent)", xlab="Horizon (years)", col="red", main="SDRs from AR(3) with exponential weights")
lines(1:N, y2, lwd=2, col="steelblue")
abline(h=0, lty=2)
abline(h=etabar1, lty=2, col="red", lwd=2)
abline(h=etabar2, lty=2, col="steelblue", lwd=2)
legend("topright", legend=c(year1, year2), col=c("red", "steelblue"), lwd=2, bg="white")
